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A PDE APPROACH TO NUMERICAL FRACTIONAL 

DIFFUSION 


RICARDO H. NOCHETTO, ENRIQUE OTAROLA, AND ABNER J. SALGADO 


Abstract. Fractional diffusion has become a fundamental tool for the model¬ 
ing of multiscale and heterogeneous phenomena. However, due to its nonlocal 
nature, its accurate numerical approximation is delicate. We survey our re¬ 
search program on the design and analysis of efficient solution techniques for 
problems involving fractional powers of elliptic operators. Starting from a lo¬ 
calization PDE result for these operators, we develop local techniques for their 
solution: a priori and a posteriori error analyses, adaptivity and multilevel 
methods. We show the flexibility of our approach by proposing and analyzing 
local solution techniques for a space-time fractional parabolic equation. 


1. Introduction 

Diffusion is the tendency of a substance to evenly spread into surrounding space, 
and is one of the most common physical processes. The classical models of diffusion 
lead to local and thoroughly studied equations. However, in recent times, it has 
become evident that many of the assumptions that lead to these models are not 
always satisfactory or not even realistic in practice. Consequently, different models 
of diffusion have been proposed, fractional diffusion being one of them. 

Capturing the essential behavior of fractional diffusion with the simplest and 
crudest models is of paramount importance in science and engineering. This al¬ 
lows for understanding of physical applications where long range or anomalous 
diffusion is considered complex phenomena in mechanics [S], biophysics |22j . 
turbulence [3T], image processing [39], nonlocal electrostatics [46], finance |50| and 
the control of devices IlISj. To understand the behavior of fractional diffusion, 
computational science is fundamental. It is one of the pillars, together with the¬ 
ory and experiments, of scientific inquiry. A carefully crafted computational model 
can replace a very expensive or unrealizable experimental setting, and it can give 
new insight into the theoretical developments of a specific discipline. The analysis 
of such computational schemes is the realm of numerical analysis, which offers a 
rigorous mathematical description of the extent to which the computer’s output 
approximates the process of interest. These crucial aspects of modern research are 
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blended together in this paper, which deals with the design and analysis of efficient 
numerical techniques for problems involving fractional diffusion. 

The mathematical structure of fractional diffusion is shared by a wide class of 
nonlocal operators |34l 154) . These operators have a strong connection with real- 
world problems, since they constitute a fundamental part of the modeling and 
simulation of complex phenomena. It is evident that the particular type of math¬ 
ematical operator appearing in applications can widely vary and that a unified 
analysis might be well beyond reach. A more modest, but nevertheless quite ambi¬ 
tious, goal is to develop computational tools and their analysis for a representative 
of a particular class. We discuss below our contributions to fractional diffusion. 

Exploiting the breakthrough by L. Caffarelli and L. Silvestre [25], we have made 
a decisive advance in numerical fractional diffusion. Although the analysis of our 
method is intricate, its implementation is done using standard components of finite 
element analysis [sniiMiiiiisniiM]- This is the main advantage of our scheme, since 
alternative approaches require less traditional techniques (such as special quadra¬ 
ture) to cope with the mathematical difficulties inherent to fractional diffusion. 
Even in ID [451166] , the integral formulation of fractional diffusion is notoriously 
difficult from the numerical standpoint due to the presence of a nonintegrable ker¬ 
nel; see m for a discussion. In contrast, our PDE approach to fractional diffusion, 
as well as the recent method of m, can handle multidimensions easily and effi¬ 
ciently - a highly desirable feature. 

Here we are interested in the design and analysis of numerical techniques to 
solve problems involving fractional powers of the Dirichlet Laplace operator (—A)®, 
s G (0,1), the so-called fractional Laplacian. Let H be a bounded Lipschitz domain 
of K" {n > 1), with boundary 917. Given a function /, we seek u such that 

(1.1) (—A)®m = / in 17. 


Our approach, however, is not particular to the fractional Laplacian and it can be 
applied to any second order, symmetric and uniformly elliptic operator [501 §7]- 
One of the main difficulties to study (1.1) is that the fractional Laplacian is a 
nonlocal operator [OH [1011031113110] ■ To localize it, Caffarelli and Silvestre showed 
in [00] that any power s < 1 of the fractional Laplacian in M" can be realized as an 
operator that maps a Dirichlet boundary condition to a Neumann-type condition 
via an extension problem on the upper half-space This result was adapted in 

[MHOIEI] to bounded domains 17, thus obtaining an extension problem posed on 
the semi-infinite cylinder C := 17 x (0, oo). This extension is the mixed boundary 
value problem: 

d'% 

(1.2) div (j/“V'^) = 0 in C, = 0 on -—= ds/on 17 x {0}, 

dv°^ 


where BlC := 917 x (0,oo) is the lateral boundary of C, and ds '■= 2^ 
a positive normalization constant that depends only on s; see [23l [25] for details. 
The parameter a is defined as a = 1 — 2s € (—1,1), and the so-called conormal 
exterior derivative of at 17 x {0} is 
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(1.4) ds{-Ayu= in rj. 


We propose, and analyze in Section the following strategy to solve 
given / we solve ([m, thus obtaining a function ; setting u : x' G i—>■ u{x') = 

{x', 0) G K, we obtain the solution of G3L The main advantage of this algorithm 
is that we are solving the local problem (1.2) instead of dealing with the nonlocal 
operator (—A)'* of (1.1). However, this comes at the expense of incorporating 
one more dimension to the problem, thus raising the question of computational 
efficiency. This motivates the use of anisotropic meshes to compensate for the 
singular behavior of '^(•, y) as y —> 0, and well as the development of a posteriori 
error estimators and multilevel methods. These are reviewed in Sections and 
respectively. To show the flexibility of our approach, in Section we consider a 
parabolic equation with fractional diffusion and fractional time derivative. 

Throughout this work 17 is an open, bounded and connected domain of K”, n > 1, 
with polyhedral boundary 917. Given the semi-infinite cylinder C = 17 x (0, oo) and 
T > 0, the truncated cylinder with base 17 and height (T is defined by := 
17 X (0,(T) with lateral boundary d^Cy := 917 x (0, (X). 

We will be dealing with objects defined in and it will be convenient to 

distinguish the extended dimension. A vector x G will be denoted by 


X = (xi,..., x„, x„+i) = (x', x„+i) = (x', y), 


with Xi G K for i = 1,..., n -I- 1, x' G K" and y G ffi. 

By a b 'we mean a < Cb, with a constant C that neither depends on a, b or 
the discretization parameters. The value of C might change at each occurrence. 


2. A PDE approach: formulation and FEM 

For a bounded domain there are several ways, not necessarily equivalent, to 
define the fractional Laplacian; see mi ISO] for a discussion. As in [50] we adopt the 
definition based on spectral theory [TO]. The operator (—A)“^ : L^(17) —>■ L^(17) is 
compact, symmetric and positive, whence its spectrum {A^ }fcGN is discrete, real, 
positive and accumulates at zero. Moreover, there exists {v^felfeeN, which is an 
orthonormal basis of L^(17) and satisfies —A(pk = XkVk in XI and = 0 on 917. 
Fractional powers of the Dirichlet Laplacian can be defined for w G C^(17) by 

OO 

( 2 . 1 ) {-AYw = '^Xlwk‘Pk, 

k^l 

where Wk = Jq w(pk- By density (—A)® can be extended to the space 

{ oo oo 

w = Wk^pk ■ y xiwi < oo 

k=l k=l 


For s G (0,1) we denote by El '^(17) the dual of IHI'’(17). 
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2.1. The Caffarelli-Silvestre extension problem. The Caffarelli-Silvestre ex¬ 
tension leads to a nonuniformly elliptic equation [53 [THl [531 [55]. We consider 
weighted Sobolev spaces with the weight \y\°‘, a S (—1,1). If C D) 

is the space of measurable functions on D such that 


GD) = [ \y\ 

Jd 


w < oo. 


m\Lmy\o 

Similarly we define the weighted Sobolev space 

H\\yr,D) = {w& L^\yr,D) : |Vn;| S L\\i 
where Vw is the distributional gradient of w. We equip H^( 


(2.3) 


\w\\mi\y\‘-,D) 




+ ll'^^lll,2(|y|“.D)) 


,D)}. 

D) with the norm 
1/2 


Since a € (—1,1) then \y\°‘ belongs to the Muckenhoupt class yl 2 (K"^^) 

IS71I73]. This, in particular, implies that H^{\y\°‘,D) with norm ( |2.3| ) is Hilbert and 
C°°{D) n H^{\y\°‘, D) is dense in H^{\y\°‘,D) (cf. [73 Proposition 2.1.2, Corollary 
2.1.6], [35] and [33 Theorem Ij). We recall the definition of Muckenhoupt classes. 

Definition 2.1 (Muckenhoupt class H2, [S7[^). Let w be a weight and TV > 1. 
We say that to G H 2 (M^) if 




(2.4) 


where the supremum is taken over all balls B in . 


< oo. 


If W e H 2 (k^), we call it an H2-weight, C 2 ,uj in \2A\ is the Al2-constant of w. 


To study problem (1.2) we define the weighted Sobolev space 
(2.5) Hl[y^,C) = {zc G H\y^,C) : zc = 0 on SlC] . 

As |60L (2.21)] shows, the following weighted Poincare inequality holds: 


( 2 . 6 ) 


\w 


U 2 (y“.C) < l|Vzz;||L 2 (y<..c), Vzz; G Hl{y°‘,C), 


then the seminorm on Hj^{y°',C) is equivalent to (2.3). For w G H^{y°',C), tinw 
denotes its trace onto H x {0} which satisfies [60l Proposition 2.5] 

(2.7) tin Hl{y°‘, C) =W{ft), || tr^ znjjHqn) < Ctr*, ||zz;l|^i 

Let us now describe the Caffarelli-Silvestre extension result for bounded domains; 
[5S1[75|- Given / G IHI“®(0), let u G ff(n) solve (1.1). If ^ G Hj^{y°‘,C) solves 
(1.2), then u = ‘^(•,0) and (1.4) holds. 

2.2. A priori error analysis. We now review the main results of m about the a 
priori error analysis of discretizations of problem (1.1). This will also serve to make 
clear the limitations of this theory, thereby justifying the quest for an a posteriori 
error analysis. In this section we assume that 


( 2 . 8 ) 


kllff 2 (n) < II A,r/zz; 11^,2(H), Vzc G i7^(H) n 


This holds if, for instance, the domain H is convex [35] - 

Since C is unbounded, problem (1.2) cannot be approximated with finite-element¬ 
like techniques. However, since the solution of problem (1.2) decays exponentially in 
y [53 Proposition 3.1], by truncating C to Cy and setting a homogeneous Dirichlet 
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condition on ?/ = fK, we only incur in an exponentially small error in terms of fK 
[501 Theorem 3.5]. If 

Hl{y°‘,Cy) := {n G 77^(y“,Cy) : n = 0 on U x {fX}} , 
then the aforementioned problem reads: find v G ^Cy) such that 


(2.9) 


y^^WvWcj) = dsif, tin (j)), 'iv & Hl{y^,Cy). 


Here (•, •) is the duality pairing between IH ®(fl) and see (2.7|. 

If and V denote the solution of (1.21 and ( |2.9[ ), respectively, then [001 Theorem 
3.5] provides the following exponential estimate 

„-vThy/4|| 


|V('^ -i;)||l 2 (j/=,c) ^ e” 


llH-=(n)- 


To study the finite element discretization of (2.9) we must understand the regu¬ 
larity of ^ and V. As [001 Theorem 2.7] reveals, the second order regularity of 
is much worse in the extended direction, namely 


(2.10) + \\dy\Ix'‘^\\L-^{y<=‘fi) < ||/||hi—( n), 

(2-11) ll‘^i/!/llL2(y/3,C) ^ ll/llL2(n)> 


with /3 > 2q;-|-1. Thus, graded meshes in the extended variable y play a fundamental 
role. Estimates (2.10)-(2.11) motivate the construction of a mesh over as follows. 
We consider a partition ly of the interval [0, fT] with points 


( 2 . 12 ) 


yk = k^M-^‘y, fc = 0,...,M, 


where 7 > 3/(1 — a) = 3/(2s). Let = {K} be a conforming and shape regular 
mesh of H, where K C K" is an element that is isoparametrically equivalent either to 
the unit cube [0,1]” or the unit simplex in K”. The collection of these triangulations 
is denoted by Tq. We construct the mesh S^y of Cy as the tensor product of 
and ly. By construction 3^y is anisotropic in the extended variable (cf. [501 160] '). 
The set of all triangulations of Cy obtained with this procedure is T. 

For G T, we define the finite element space 

(2.13) Y{sry) = {W G C°{Cy) : W\t G Vi{K) ® Pi(/) VT G Wjr^ = 0} , 

where T d = dyCy U H x {fT} is the Dirichlet boundary. The set Vi{K) is either 
the space Pi (AT) of polynomials of total degree at most 1, when the base K of 
an element T = K x I is simplicial, or the space of polynomials Qi(Ar) of degree 
not larger than 1 in each variable provided K is an n-rectangle. We also define 
U(.%) = traV(c^). Note that #^y = Mand that « M” implies 

« M”+b 

The Galerkin approximation of (|2.9|) is the function V 57 G 'Vi^y) such that 


(2.14) f y“VE 9 -^VW = 4(/,trnW), VW G V(5^y). 

Jc^ 

Existence and uniqueness of t /37 immediately follows from V(,i5^) C Hy{y°‘,Cy) and 
the Lax-Milgram Lemma. It is trivial also to obtain a best approximation result d 
la Cea. This reduces the numerical analysis of (2.14) to a question in approxima¬ 
tion theory which in turn can be answered by the study of piecewise polynomial 
interpolation in Muckenhoupt weighted Sobolev spaces; see [001103]. Exploiting the 
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Cartesian structure of the mesh we can handle the anisotropy, construct a quasi 
interpolant 11^^. : L^(Cy) —)■ V(^), thus extending [35], and obtain 


||z) - < hK\\'^x'V\\L^(yO.^ST) + hl\\dyV\\L2(^ya. ^St) ^ 

\\dxjiv - J\^^v)\\L2l^yc^T) < hxWVx'dxjV\\L2l^y,^^ST) + hl\\dydxjV\\L2(^ya,^ST)^ 


with j = 1,... ,n + 1; see [601 Theorems 4.6-4.8] and [63]. However, since (2.11) 
implies ^ ,Cy) the second estimate is not meaningful for j = n + 1. We 

must measure the regularity of ^yy with a stronger weight and compensate with 
a graded mesh. This makes anisotropic estimates essential. These considerations 
allow us to obtain l [5ni Theorem 5.4] and [501 Corollary 7.11]): 


Theorem 2.2 (a priori error estimate). Let SLy G T and V(.^) be defined by 
(2.13). IfVs^ GY{t7y) solves (2.14), then we have 




-l/(n+l)| 




where y « \og{fi^3y). Alternatively, if u denotes the solution of (1.1), then 

lk-troH5yllH=(n) < |log(#^r)r(#^r)-'^^”+'^ll/llHi-=(0). 


Remark 2.3 (domain and data regularity). The results of Theorem 2.2 hold only 
if / G IHI^“®(H) and the domain H is such that (2.8) holds. 


Remark 2.4 (quasi-uniform meshes). Let ^ solve (1.2), and V be the solution 
of (2.14), constructed over a quasi-uniform mesh of size hg. If / G IHI^“®(H) and 
y ~ I \oghg\, then for all e > 0 


(2.15) \\^{^ -Vg,)\W(ygc,)<h^i^\\ 

where the hidden constant blows up if e tends to 0. 


qn), 


Remark 2.5 (Case s = |). If s = then there is no weight in (2.15) and we 
obtain the optimal estimate 


l|V('^-y^,)llz,qc,)<h^||/||^V.(^). 

2.3. Numerical illustrations. Let us present several numerical examples. These 
have been implemented with the deal. II library mm- Integrals are evaluated 
with Gaussian quadratures of sufficiently high order and linear systems are solved 
using CG with ILU preconditioner and the exit criterion being that the f^-norm of 
the residual is less than 10“^^. 


2.3.1. Quasi-uniform meshes. According to m, dy^f/ ~ y “, which formally im¬ 
plies dy'^ G ,C) provided / G IHI^“'*(r2). In this sense (2.15) seems sharp 

with respect to regularity even though it does not exhibit the optimal rate. We 
verify this with a simple numerical illustration in Figure [2T] 

Let n = 1, H = (0,1), s = 0.2, / = sin(7ra:), then u{x) = sin(7rx), and the 
solution 'f/ to (1.2) is 'f/{x,y) = 2^“®7r®r(s)“^ sin(7ra:)Ars(7rj/), where by we 
denote the modified Bessel function of the second kind; [551 §2.4]. Figure [2T] shows 
the rate of convergence for the i7^(2/“,C9^)-seminorm. Estimate ( |2.15| ) predicts a 
rate of . In other words the rate of convergence, when measured in terms of 

degrees of freedom, is {ff^y)~^'^~^ 


which is what Figure 2.1 displays. 
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Figure 2.1. Rate of convergence for quasi¬ 

uniform meshes: s = 0.2, n = 1. 


m,n gN. 


2.3.2. Graded meshes: a square domain. Let 12 = (0,1)^. Then 

ipm,n{xi,X2) = sin(TO7ra:i) sin(n7ra;2), (m^ -|- n^) , 

If/(aii,a: 2 ) = (27r^)® sin(7ra:i) sin( 7 rai 2)5 then u(a:i,a; 2 ) = sin( 7 ra;i) sin( 7 ra; 2 ), by (2.1) 
and ^{xi,X 2 ,y) = 2 ^“®/^ 7 r®r(s)“^ sin( 7 ra:i) sin( 7 ra: 2 )y*Ris(v^ 7 r 2 /) [201 (2.24)]. 

We construct a sequence of meshes {■Shy^}k>i, where =% is obtained by uniform 
refinement and is given by (2.12) with parameter 7 > 3/(1 — a). On the basis 
of Theorem 2.2 the truncation parameter % is chosen by 

% > ^(logC-log(#5^y,_J-i/3). 

With this type of meshes, 

||u - trn <11^- ^ I log(#^rJr • (#^7.)”'/', 

which is near-optimal in but suboptimal in u, since we should expect (see m) 
Ik - tro F^,.J|h 70 ) < < (#17^J-(2-«)/3. 

Figure [ 2 . 3 .2 1 shows the rates of convergence for s = 0.2 and s = 0.8 respectively. 


In both cases, we obtain the rate given by Theorem 2.2 




Degrees of Freedom (DOFs) Degrees of Freedom (DOFs) 


Figure 2.2. Computational rate of convergence for a square and 
graded meshes. The left panel shows the rate for s = 0.2 and the 
right one for s = 0.8. The experiment al ra te of convergence is 
in agreement with Theorem 


2.2 
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Figure 2.3. Computational rate of convergence for a circle with 
graded meshes. The left panel shows the rate for s = 0.3 and the 
right one for s = 0.7. The experiment al ra te of convergence is 
j agreement with Theorem 


2.2 


2.3.3. Graded meshes: a circular domain. If = {|a;'| G : \x'\ < 1}, then 

(2.16) (pjji .fi(rj 0^ — Jm(^jm,nX^ (^m,n COs(777-0) T Bm.n sin(77I^)) , Xm.n — jm,m 

where Jm is the m-th Bessel function of the first kind; jm,n is the n-th zero of J„ 
and Am,n, Bm,n are normalization constants to ensure ||i^m,n||L 2 (a) = 1- 


If / = then (2.1) and [501 (2-24)] show that u = v?i,i and 

'^(.r,e,y) = 2^~’'T{s)~^{Xi^iy/^(pi^i{r,e)y''Ks{V2Try). 

We construct a sequence of meshes {Ay^}k>i as in p.3.2 With these meshes 


(2.17) 




- 1/3 


which is near-optimal. Figure 2.3 shows the errors of for 


s = 0.3 and s = 0.7. The results, again, are in agreement with Theorem 2.2 


3. A posteriori error analysis AND ADAPTIVITY 

Starting with the pioneering works [SIT] in the late 1970’s, adaptive finite ele¬ 
ment methods (AFEM) have been the subject of intense research. This is because 
AFEM yields optimal performance in situations where classical FEM cannot. In our 
problem we are incorporating one extra dimension and so AFEM is essential to re¬ 
tain computational efficiency. In addition, the a priori theory requires / G IHI^“®(f2) 


and (2.8). If either of these does not hold ^ may have singularities in the x'- 


variables and exhibit fractional regularity. This would not allow us to attain the 
almost optimal rate of convergence of Theorem [2.2| a quasi-uniform refinement of 
n would not result in an efficient solution technique. An adaptive loop driven by an 
a posteriori error estimator is fundamental to recover optimal rates of convergence. 
Let us explore this now and develop a new estimator. 

3.1. Residual estimators. A residual error estimator uses the strong form of the 
local residual to measure the error. Let T G S^<y and v be the unit outer normal to 
T. Integration by parts yields 


• VIE = 


/ IEy“VF7^-u - / div(y“VF5y)IE. 

IdT JT 
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Since a S (—1,1) the boundary integral is meaningless for y = 0. Even the very 
first step in the derivation of a residual a posteriori error estimator fails! There is 
nothing left to do but to consider a different type of estimator. 


3.2. Local problems on stars over isotropic refinements. Following [51 [5S] we 

can construct, over shape regular meshes, an error estimator based on the solution of 
small problems on stars. Its construction and analysis is similar to the developments 
of §3.3[ so we skip details. Under some assumptions this estimator is equivalent 
to the error up to data oscillation; see [62] ■ We designed an adaptive algorithm 
driven by this error estimator on shape regular meshes [221 [56], and we illustrate 
its performance with a simple but revealing numerical example. 

Let n = (0,1) and s e (0,1). If f{x') = sin(7ra;'), then u{x') = sin(7ra;'), and 

‘^{x',y) = 2^“®7r®r(s)“^ sin(7rx')i4rs(T?/) solves (1.2|. We point out that we are 
solving a two dimensional problem so the optimal rate we expect is 0(# 


Figure |3.1| shows the experimental rate of convergence for s = 0.2 and s = 0.6 
which, as we see, is and coincides with the suboptimal one obtained 

with quasi-uniform refinement [601 §5.1]. These numerical experiments show that 
adaptive isotropic refinement cannot be optimal, thus justifying the need to intro¬ 
duce cylindrical stars together with a new anisotropic error estimator, which will 
treat the ^'-coordinates and the extended direction y separately. 



Figure 3.1. Computational rate of convergence #(.3^) for an 
isotropic adaptive algorithm for n = 1, s = 0.2 and s = 0.6. 


3.3. Local problems over cylindrical stars. Inspired by [6] [56], we deal with 
the anisotropy of the mesh in the extended variable y and the coefficient j/“ by con¬ 
sidering local problems on cylindrical stars. The solutions of these local problems 
allow us to define an anisotropic a posteriori error estimator which, under certain 
assumptions, is equivalent to the error up to data oscillation terms. 

Given a node z on the mesh we exploit the tensor product structure of 
and we write z = (z', z") where z' and z" are nodes on the meshes and ly 
respectively. For K G we denote by C^i{K) and t^{K) the set of nodes and 
interior nodes of respectively. We set 

KeSn KeSn 

The star around z' is Sz> = Uicgz' K C fl, and the cylindrical star around z' is 
Cz' --[JiT G -.T = K X I, K^z'} = Sz' x (0, fT) C Cy. 
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For each node z' S ^A£(^) we define the discrete space 

W(C^/) = {Wg C{Cy) : W\t G V 2 {K) 0 P 2 (/) 'iT = K x I 
^bc^/\ax{o} = 0} j 

where, if if is a quadrilateral, V 2 {K) = Q 2 (if) — the space of polynomials of degree 
not larger than 2 in each variable. If if is a simplex, V 2 {K) = P 2 (if) 0 ]B(if) where 
P 2 (if) is the space of polynomials of total degree at most 2, and IB(if) is the space 
spanned by a local cubic bubble function. 

For each cylindrical star C^/ we define 7]^' G W{Cz') to be the solution of 

^ i/“V7;,-VVF = d,(/,trnVF)-^ VIF G W(C,0- 

We finally define the local indicators and global error estimators S'sTy as follows: 

(3.1) = <^’ 4 = E 

2'63v;(5h) 

We can obtain a local lower bound for the error without any oscillation term and 
free of any constant. 

Theorem 3.1 (localized lower bound). Let v G ,Cy) and V G solve 

(2.9) and (2.14) respectively. Then, for any z' G ‘JifSTvi), we have 

For z' G iV((.l5n), we define the local data oscillation as 

(3.2) oscAff ■■= dshlmf - UIk ■■=£/, 

whence the global data oscillation reads 

■= E osc^'(/)^- 

z'e^iSa) 

To mark elements for refinement we use the total error indicator 

(3.3) T^AV^„SA~{'^i+oscAfff'^ Vz'GAiAi). 

Let = {Szi : z' G for any C set 

1/2 

(3.4) I E 

Under certain assumptions we can bound the error by the estimator, up to 
oscillation terms; see [30] for details. 

Theorem 3.2 (global upper bound). Let v G H} {Cy,yA and G N{SLy) solve 
( [^ and ( |2l4t > respectively. The total error estimator T 3 ri^{V^y,,,JAh), defined in 

(3.4) satisfies 

||V(u - VyrA\L^y-‘,Cr) ^ rSniVsKy, AyA- 
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3.4. Numerical experiments. We now illustrate the performance of the a pos¬ 


teriori error estimator (3.1 1 . We use an almost standard adaptive loop 
(3.5) 


SOLVE ^ ESTIMATE ^ MARK ^ REFINE. 


where the modules in (3.5) are as follows: 


SOLVE: Given a mesh S/'y we compute the solution of (2.14). 


ESTIMATE: Given we calculate the local error indicators (3.1 1 and the local 


to construct the total error indicator (^3.3 


oscillations 

= ESTIMATE(V5-,,5'^). 

MARK: Using Dorfler marking [53] with parameter 0 < 6 * < 1 we select a set 
of minimal cardinality that satisfies t^ 

REFINE: We generate a new mesh by bisecting all the elements K G 
contained in ^ based on newest-vertex bisection method; jSSlIMj- We choose 
the truncation parameter as (T = l-F | log(#.^) [601 Remark 5.5]. We set M « 


and construct ly by the rule (2.12). The new mesh = REFINE(.^) 
is obtained as the tensor product of and Xy. 

3.4.1. Smooth but incompatible data. The example of |6()1 §6.3] shows that Theo- 
is sharp: / G is necessary to obtain an optimal rate of convergence 


rem 


2.2 


with a quasiuniform mesh in the x'-direction. A certain compatibility between the 
data and the boundary condition is necessary. Moreover, [Sni §6.3] shows how, in 
some simple cases, one can guess the singularity and a priori design a mesh that 
captures it and recover the optimal rate of convergence. This is not always possible 


and here we show that the estimator (3.1) automatically produces a sequence of 
meshes that yield the optimal rate of convergence. Consider U = (0,1)^ and / = 1. 
From (2.2) we see that functions in IHI^“®(r2), for 1 — s > must have a vanish¬ 
ing trace. Therefore, 


if s < i, / ^ 


irl-s 


(11) and Theorem 


2.2 


cannot be invoked. 


Nevertheless, as Figure 3.2 shows, we recover the optimal rate. 



10^ 10'' 
Degrees of Freedom (DOFs) 



10^ 10'* 
Degrees of Freedom (DOFs) 


Figure 3.2. Computational rate of convergence of AFEM for 
§ [3AT| n = 2 and s = 0.2, 0.4, 0.6 and s = 0.8. The left 
panel shows the the error vs. the number of degrees of freedom, 
the right one the total error indicator. We recover the optimal 
rate #(f5y)“^/^. For s < 5 , the right hand side / = 1 ^ IHI^“®(S1) 
and a quasiuniform mesh in U does not deliver the optimal rate of 
convergence [601 §6.3]. 
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3.4.2. L-shaped domain with incompatible data. We now combine the singularity 
introduced by the data incompatibility of § |3.4.1| and the effect of a reentrant 
corner. Consider = (—1,1)^ \ (0,1) x (—1, 0) and / = 1. As Figure 3.3 displays, 
we recover the optimal rate of convergence for all possible cases of s. 



10® lO" 

Degrees of Freedom (DOFs) 



Degrees of Freedom (DOFs) 


Figure 3.3. Computational rate of convergence of AFEM on the 
smooth but incompatible right hand side over an L-shaped domain 
of § |3.4.2| for n = 2 and s = 0.2, 0.4, 0.6 and s = 0.8. The left panel 
shows the error vs. the number of degrees of freedom, whereas the 
right one that for the total error indicator. In all cases we recover 
the optimal rate 


4. Multilevel methods 


Since we increase the dimension by one in the approximation of the nonlocal 


operator (—A)® by the local problem (1.2), it becomes essential to develop efficient 


linear algebra methods. It is known that multilevel methods are the most efficient 
techniques for the solution of discretizations of PDFs [191 [20l|43l[75]. Multigrid 


methods for equations of the type (1.2), however, have not been explored. 


Using the multilevel framework of mini [75], the Xu-Zikatanov identity m\ 
and exploiting the fact that |y|“ G A 2 (M."~^^), which turns out to be critical, we 
derived in [29] an almost uniform convergent multilevel method to solve (2.14). As 


we have shown the mesh in the extended dimension must be graded towards the 
bottom of the cylinder thus becoming anisotropic. We apply line smoothers over 
vertical lines in the extended domain and prove that the corresponding multigrid 
V-cycle converges nearly uniformly, i.e., the contraction factor depends linearly on 
the number of levels, and thus logarithmically on the problem size. 


4.1. Multilevel decomposition and multigrid algorithm. We follow [13] [14] 
to present a multilevel decomposition of V(c^). To simplify the analysis and imple¬ 
mentation, we consider a sequence of nested discretizations constructed as follows: 
We introduce a sequence of nested uniform partitions of the unit interval {Ifc}, with 
mesh points pi^k, for I = 0,, Mk and fc = 0,..., J. For 7 > 3/(1 — a) the family 
{Ik} is given by 


(4.1) yi,k = ‘yylk, l = 0,...,Mk. 

For k = 0,..., J, is obtained by uniform refinement. The mesh is the 
tensor product of A{i^k and Ik, given by (4.1). 
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If Vfc := we have the sequence Vq C Vi C • • • C Vj = V, and 

a space macro-decomposition V = introduce a space micro¬ 

decomposition: For j = 1,..., ?V4 let Ifc be a subset of {1,2,..., !A4'}, and assume 
that ={1,2,..., !A4'}. The sets may not be disjoint but the cardinality 

of their intersection is bounded independently of J and IA6 c- If the nodal basis of Vfc 
is 4>k,i, we define 'Vkj = span{(/)fc ^ : i G ^kj} and we have the micro-decomposition 


(4.2) 


.7 314 


j' 


k—Oj—1 


With this notation we define a symmetric V-cycle multigrid method as in na 
Algorithm 3.1], with m > 1 pre and post smoothing steps. When tti = 1, it is 
equivalent to the application of successive subspace corrections to the decomposition 


(4.2) with exact sub-solvers so that the V-cycle multigrid method has a smoother at 


each level of block Gauss-Seidel type nans]- In particular, the nodal decomposition 
Ik,j = {j} yields a point-wise Gauss-Seidel smoother. If the indices in Ikj are such 
that the corresponding vertices lie on a vertical line, we obtain line smoothers, 
which are essential to handle anisotropy. 


4.2. Analysis of the multigrid method. The success of multigrid methods for 
uniformly elliptic operators is because smoothers are effective in reducing the high 
frequency components of the error while coarse grid corrections reduce the low 
frequency ones. However, their effectiveness depends on several factors such as 
the shape of the mesh. A key ingredient in the design and analysis of a multigrid 
method on anisotropic meshes is the use of line smoothers n El min]- 

When solving (2.14) on graded meshes, the approximation on the coarse grid 


is dominated by the larger mesh size in the cc'-direction and thus the coarse grid 
correction cannot capture the smaller scale in the y-direction. This is why line 
smoothing, i.e., solving sub-problems restricted to one vertical line is necessary. 
Owing to the nature of the decomposition, the smoother requires the evaluation of 
the inverse of the operator over a vertical line. This can be efficiently realized since 
the corresponding matrix is tridiagonal. 

Our analysis hinges on IZ21, which requires stability of the micro-decomposition. 

Lemma 4.1 (nodal stability and anisotropic inverse inequalities). Let G 
with quasiuniform and Xy graded so that (4.1) holds. If v = O ^ 

then 


(4.3) 


Wj 

E"^ 

i=i 




5V7j 

- E 

7 = 1 


OllL2(y=,CH' 


Moreover, for T = K x I G Sly we have the following local inverse inequalities 
(4.4) \\\7,c'V\\L^(^ya.^T) ^ hj/\\v\\L 2 (^yo.^y^, ||9y'y||L2(y°.T) ^ ^ ^\\v \\^2 (^yo. 

We now examine the V-cycle multigrid method applied to the decomposition 


(4.2) with exact sub-solvers on Vfej, i.e., with line smoothers; see [1^1 §111.12] 
and [74]. A key observation in favor of subspaces follows. 
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Lemma 4.2 (nodal stability of y-derivatives). In the setting of Lemma 4-.1 we have 

9v{j 9v[j 

(4.5) ^ \\dyVj\\li(yc^c^) < ^ \\dyvf'^ 


i=i 


i=i 


We use the quasi-interpolant of [63], Lemmas 4.1 and 4.2 and obtain nearly 
uniform convergence of the V-cycle multigrid method. We follow laizei. 


Theorem 4.3 (convergence of multigrid with line smoothers). The symmetric V- 
cycle multigrid method with line smoothing converges with a contraction rate 


5<\- 


1 

l + CJ’ 


where C is independent of 9fj and depends on y°‘ only through ^ 2 ,^°. 


4.3. Numerical illustrations. We consider two cases: 

• n = 1, n = (0,1), u = sin(37ra;), 

• n = 2, n = (0,1)^, M = sin(27ra:i) sin(27ra;2), 

and 9^ = 1, which has been chosen, as discussed in |60j . to capture the exponen¬ 
tial decay of the solution. All of our algorithms are implemented based on the 
MATLAB® software package i FEM [55] • 

Table |LT| shows the number of iterations for the one and two dimensional prob¬ 
lems, respectively. As we see, the method converges almost uniformly with respect 
to the number of degrees of freedom. 



DOFs 

s = 0.15 

s = 0.3 

s = 0.6 

s = 0.8 

16 

289 

7 

6 

5 

5 

i 

32 

1,089 

9 

9 

6 

6 

r 

64 

4,225 

10 

10 

6 

6 

I 

128 

16,641 

10 

11 

6 

6 

1 

256 

66,049 

11 

10 

6 

6 

1 

512 

263,169 

11 

10 

6 

7 


hSa 

DOFs 

s = 0.15 

II 

O 

CO 

s = 0.6 

00 

o 

II 

CO 

1 

16 

4,913 

8 

7 

6 

5 

1 

32 

35,937 

11 

8 

6 

6 

1 

64 

274,625 

12 

9 

6 

6 

I 

128 

2,146,689 

13 

9 

6 

6 


Table 4.1. Number of iterations for a multigrid method for 
{—Ayu = f with n = 1 (top) and n = 2 (bottom) using a line 
smoother in the extended direction. The mesh in LI is uniform 
with mesh size h^^, whereas that in the extended direction (0,fF) 
is graded according to (4.1). 


5. Space-time fractional diffusion problems 

We now review the numerical approximation of an initial boundary value problem 
for a space-time fractional parabolic equation developed in |S5j. Given s € (0,1), 
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7 G ( 0 , 1 ], a function /, and an initial datum uq, we seek u such that 
(5.1) d^u + {—A)^u = f inilx( 0 ,T), u{0) = uq in 


If 7 G (0,1), df denotes the left Caputo derivative of order 7 , defined by 


(5.2) 


dfu{x,t) := 


du{x, r) 


r(l-7)io dr 


dr, 


where F is the Gamma function. For 7 = 1 , we consider the usual derivative 9*. 
Fractional derivatives and integrals, like the Caputo derivative, are a powerful 


tool to describe memory and hereditary properties of materials mi- Problem ( |5.1[ ) 
is motivated by anomalous diffusion processes (subdiffusion) and models dynami¬ 
cal systems with chaotic motion |BS], fractional PID controllers [5B], subdiffusion 
phenomena in highly heterogeneous aquifers [55] and plasma turbulence |32| . 


Problem (5.1 1 is complicated due to the nonlocality of the fractional time deriv¬ 


ative jUKini besides the presence of the nonlocal space operator. We overcome the 
latter with the Caffarelli-Silvestre extension of § 2.1 and rewrite (5.1) as an elliptic 
problem with dynamic boundary condition: 


(5.3) 


div (i/“V'^) = 0 in C X (0, T), ^ = uq, on fl x {0}, t = 0 

= 0 on 9z,C X (0, T) dsdffif + = dsf on {Vt x {0}) x (0, T). 


There are several approaches via finite differences, finite elements and spectral 
methods to treat the Caputo derivative of order 7 . We refer to m § 1 ] for an 
overview. In |521153j a finite difference scheme is proposed, which has a consistency 
error where r denotes the time step. This estimate, however, requires 

a rather strong regularity assumption in time |55[ 162] . In |62j , we examined the 
behavior of dtu and duu when t 0 and derived realistic time-regularity estimates 
for u; see also [SaEHl- With these refined results we analyzed the truncation error 
and showed discrete stability. The latter leads to an energy estimate. 


5.1. The Caffarelli-Silvestre extension. We define 

W := {w G L^{0,T;L'^{n)) n : dfw G T; H-"(0))}, 

V := {zc G L\0,T-Hi{y‘^,C)) : Of tmw G L\0,T;U-%n))}. 


Thus, given / G L^(0, T; IHI“'’(n)), a function u GW solves (5.1) if and only if the 
a-harmonic extension G V solves (5.3). A weak formulation of (5.3) reads: Find 
G V such that tro (0) = uq and, for a.e. t G (0, T), 

(5.4) (tra97'^,tro(/))-h /'y“V'^V(/)= (/,trn^) \/cj) G 

as Jc 


Remark 5.1 (dynamic boundary condition). Problem (5.4) is an elliptic problem 
with a dynamic boundary condition: = dgf — dsirndf^ on x {0}. Its 

analysis is slightly different from the standard theory for parabolic equations. 

Given s G (0,1), 7 G (0,1], a forcing term / G L^(0, T; IHI“®(fl)) and Uq G T^(II), 
problems (5.1) and (5.3) have a unique solution; see [511 Theorem 2.6]. 
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5.2. Time regularity. As in the elliptic case, in what follows we assume that (2.8) 
holds. We refer to [55] for a complete space regularity analysis of problem (5.1). 
We focus here on time regularity. 

For 7 = 1, we demand sufficient time regularity of the right-hand side together 
with compatibility conditions for the initial datum uq. We express this as 

(5.5) tvadu^ eL\0,T;M-^{n)). 


For 7 S (0,1), (5.5) is inconsistent with the solution of (5.1). Properties of the 
Mittag-Leffler function show that (5.5) never holds if uq ^ 0. Time derivatives of 
u are unbounded as t 0. In particular, duu{x',t) ^ L^(0, T; ]HI“'^(f2)). However, 


/o+ 


r||5ttM(-,t)||^_,(f2) dt 


is finite provided cr > 3 — 27 . For this reason, when 7 G (0,1), we assume 

G L"( 0 ,T;H-"(H)) a > 3 - 27, 

which is a valid assumption provided A{uq, /) < 00 , where 

( 5 . 6 ) A{uo,f) = ||uo||H'>(n) + ll/llff2(o,T;H-7n))- 

Theorem 5.2 (time regularity: 7 G (0,1)). If f G i?^(0, T; ]HI“'*(il)) and uq G 
IHI®(H), then for t G (0,T], the solution of (5.1) satisfies 

\\dtu{-,t) - ^^■u(-,t)||H-»(n) < t'^~^A{uoJ), 

where i5^u(-, t) = t~^ {u{-,t) ~ u{-, 0)). Moreover, 

\\t'^^'^dttu\\ L2(0.T;H-7n)) ^A{uo,f), 

where cr > 3 — 27. The hidden constant is independent of t but blows up as ^ f 0 . 

5.3. Time discretization. Let /C G N denote the number of time steps. The 
time step is r = T/JC, and, for 0 < k < 1C, tk = kr and Ik = {tk,tk+i]- For 
4> G C{\f),T],X) we denote (j)’^ = 0(tfc) and (j)'^ = Moreover, 




k\\2 




For W'^ C A' we define, for A: = 0,..., /C — 1, = r — W^). 


5.3.1. Time discretization for 7 = 1. We apply the backward Euler scheme to (5.4) 
for 7 = 1: determine W = C such that 

(5.7) trn = uo, 

and, if f^~^^ = /(t^+^) for fc = 0,...,/C — 1, then 1/^+^ G H\{y°‘,C) solves 

(<5Urnr'=+\troW)z,2(a) + J- [ = (7+\tro W), VW G 7?i(y“,C). 

as Jc 

Define by = tra G IHI®(D), which is a piecewise constant (in time) ap¬ 
proximation of u, solution to problem (5.1). Note that (|5.7|) does not require an 


extension of uq. The stability of this scheme is elementary |62[ Lemma 3.3]: 

II trn l^’^||foo(L2(n)) + ll^"'ll^2(^i (^c^c)) ^ II'“oIIl2(q) -|- \\r\\%(m-‘’(n))- 
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5.3.2. Time discretization for 7 G (0,1). We now discretize df for 7 G (0,1). We 
consider the scheme proposed in [521 [53] but resort to the regularity results of 
Theorem 5.2 Using (5.2) and the Taylor formula with integral remainder yields 

I ^ 

(5.8) dfu{-,tk+i) = 


1 


r(2-7) 


k 

u{-,tk+l-j) — u{-,tk-j) I 

3=0 


for 0 < fc < /C — 1 , where 


(5.9) Oj = {j + ly -j^ 




g {tk+1 — t)' 


:R{-,t) dt, 


and 

(5.10) 


R{-,t) = dtu{-,t) - {u{-,tj+x) - u{-,tj)) yt G Ij. 


Notice that from (5.9) we deduce that aj > 0 for all j > 0 and 1 = oq > ai > 02 > 

• • • > ttj, limj_>.oo aj =0. 

5.3.3. Consistency. We estimate rl) using a cancellation property: R, defined in 
(5.10), has vanishing mean in Ij whence we can write 


^_ y 

^ r(l-7)^^. 


i'lp.yit) dt, 


with 'ip.y{t) = (t/c+i — t) ^ and fp-yit) dt. Applying [521 Lemma 2.1] yields 

lll'glL2(o.T;H-7n)) ^ IIV'7 “''^gUdO,r)|l-R^||L2(o,T;H-'>(n))) 

which reduces the estimation of the residual to deriving suitable bounds for each 
term on the right hand side of this expression. The regularity results of Theorem 
combined with the special structure of the kernel ipry, yield an estimate for rl). 


5.2 


This estimate, although giving lower rates of convergence than [521 (3.4)], takes into 
account the correct behavior of the solution and the singularity of its derivatives as 
110 [521 Proposition 3.6]: 

l|l'7llL2(o.T;H-7n)) ^ A{uo, f) 0 < 9 < 1/2. 

The hidden constant is independent of the data and r but blows up as 0 f g 


5.3.4. Stability and energy estimates. We apply (5.8). If 7 G (0,1) and (jC C L^(fl) 
we define the discrete fractional derivative, for A: = 0,..., /C — 1, by 


r (2 - := y ^<51 


3=0 


/.k+l-j _ 


^fe +1 

T'y 


k-1 


-E 

3=0 


Qj °J+1 j.k-3 _ f±^0 




T'1 


The implicit semi-discrete scheme to solve (5.1) reads: Set U° as in (5.7) and, for 
fc = 0,...,/C — 1, G H\{y°‘,C) satisfies, for every W G H\{y°‘,C), 

(5.11) {5^ trn trn W)^^) + j- ^ y^VV’^+^VW = {y\W). 

Let I'^w be the Riemann-Liouville integral of order a of w evaluated at t = T. 
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Theorem 5.3 (stability for 7 S (0,1)). The implicit semi-discrete scheme (5.11) 
is unconditionally stable and satisfies 


(5.12) /'-^||trnTnii.(a) + l|V^' 


t||2 




< ’^I|mo||l2(q) + ||/"’||£2(]H[-.(n). 


Deducing an energy estimate for problem (5.1) is cumbersome due to the nonlo¬ 
cality of the fractional time derivative. A key ingredient in deriving such a result 
is an integration by parts formula, which for a function not vanishing at t = 0 
and t = T involves boundary terms that need to be estimated; for a step in this 
direction see Eziini]. The energy estimate ( |5.12[ ) yields: 

Corollary 5.4 (energy estimate for u). Let 7 G (0,1). Then, 

(5.13) "’'ll“llL 2 (n) + ll^llL 2 (o_T;H=(n)) — ^II^oIIl 2 (q) -I- ||/|li 2 (o_T;H- 7 n))' 

Remark 5.5 (limiting case). Given g G LP(0,T), we have g —>■ g in LP{0,T) as 
CT I 0 ; see [701 Theorem 2.6]. This implies that, taking the limit as 7 f 1 in ( jAlsl ), 
we recover the well known stability result for a parabolic equation, i.e.. 


(5.14) 


li“(o.r;L2(a)) + l|n||i2(o^T;H'’(n)) ^ 


2 

L 2 (n) 


lz,2(o,T;H- = (n))' 


This allows us to unify the estimate of Corollary 5.4 for all 7 G (0, Ij. 


5.4. Truncation. We now study the space discretization of (5.4). Since ^(t) 
decays exponentially in the extended direction y, for a.e. t G (0,T), we truncate C 
to Cy for a suitable T and seek solutions in Cy. Define, for 7 G (0,1], 

(5-15) K^{uq, f ) •.= "''ll'*^o|lL 2 (n) + ll/llL 2 (o,T;H-'’(n))’ 

where, by Remark [s^ 1° is the identity. For 7 G (0,1] and s G (0,1), we have 

(5.16) IIV'^|U2(o,T;L=to<..ax(r.oo))) < e-^^/2A^(Mo,/), 


where fk > 1 and denotes the solution to (5.4). 

As a consequence of (5.16), we can consider the truncated problem 

(5.17) 


div (j/“Vw) = 0 in Cy X (0, T), u = 0 on {dyCy U Lty) x (0, T) 
dsdf trn v dfv = dsf on (fl x { 0 }) x ( 0 , T), 


with the initial condition u = mq on D x {0} and f = 0, and where Lly = Ll x {T} 
with fk > 1 sufficiently large. Upon defining 

Vy = [w& L\Q,T-Hl{y^,Cy)) : 97 ti'o £ L^0,T-,U-^{n))} , 

we understand ( |5.17 ) as: seek v GVy such that tr^ u(0) = Uq and for a.e. t G (0, T), 

(5.18) {df tin V, trn (j>)^ f y°^WvV(j) = if, trn (f), Vcj) G Hl{y°‘,Cy). 

S J Cy 

If solves (5.4), V solves (5.18) and A-y{uo,f) is given in ( |5.15 1, then we have 
the following exponential convergence result for every 7 G (0,1] and fk > 1: 

(5.19) /i-T'll tro('^ - v)\\l2^n) + II“ ^^)IIl 2 (o,T;LR y“.Cy)) (mq, /). 
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5.5. Fully discrete scheme. We describe a fully discrete scheme to solve ( |5.18 ). 
The space discretization hinges on the FEM discussed in Sectio n T he discretiza¬ 
tion in time uses the schemes proposed in 15.3.1 for 7 = 1 and in j ]5.3.2 for 7 S (0,1). 

The fully discrete scheme computes the sequence C an approxima¬ 

tion of the solution to (5.18) at each time step. We first initialize the scheme: 

(5.20) 

where = Go 'Ka-, 'H.a. denotes the a-harmonic extension onto Cy and Ggy 
the weighted elliptic projection of [H2 §4.3]; notice that tro = iTQ,GgyV[Q). 
For fc = 0,... ,/C — 1, let G satisfy for all W € 

(5.21) (S'^tTnV^+\tTnW)mn) + ^J^ = {f+\trnW). 

An approximate solution to problem (5.1 1 is given by = tr^ Vg^- 

The discrete scheme ( |5.20 )-(5.21 1 is unconditionally stable for all 7 G (0,1]: 

where 1° is the identity according to Remark 5.5 (case 7 = 1 ). 

We define B{uo,f) := |luo||Hi+ 3 '>(a) + ll/|t=o||Hi+ 7 n) + ll/lliu^(o,T;Hi-(i- 2 A‘) 7 n)) 
for /r > 0. The error estimates for (5.20)-(5.21) read as follows. 

Theorem 5.6 (error estimates for 7 G (0,1)). Let 7 G (0,1), v and solve 
(5.18) and (5.20)-(5.21), respectively. If A(uo, f), B{uo, f) < 00 and SLy satisfies 
(2.12), then 

trn(u" - R5,)lli=(n)]^ < r^A{uoJ) + \ logiVp^lV^ 6 (uo,/), 


and 


y - y^WtHHUygC,)) < 


\\ogN\^N^BiuoJ), 


where Q < 9 < I and A is defined in (|5.6[). 


We finally observe that exploiting the regularity (5.5), we can derive optimal 
(linear) error estimates in time for 7 = 1 . 
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